{"cells":[{"cell_type":"markdown","source":"# Lecture 12: Model evaluation and comparison (3) & GLM\n\n## Instructor： 胡传鹏（博士）[Dr. Hu Chuan-Peng]\n\n### 南京师范大学心理学院[School of Psychology, Nanjing Normal University]","metadata":{"id":"9790154EC7D9411EB05D6EE5D11C19FC","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##  Recap: Model evaluation & comparison (3)\n\n### 模型评估/比较与选择(Model evaluation/comparison & selection)\n\n- 过拟合\n\n- 欠拟合\n\n","metadata":{"id":"460268A615ED4EDF8FD19AED876C3931","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"### 常用的模型评估与比较的统计指标可被分为三类。 \n\n- 模型的拟合优度(Goodness of fit)\n\n-- 拟合优度(Goodness of fit)\n\n-- 均方差(Mean squared error, MSE)\n\n-- 对数似然(Log-likelihood)\n\n\n","metadata":{"id":"908B04663AD94FAE8A94522862DA6885","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"* 对新数据的预测准确性\n\t* 交叉验证\n\n* 信息准则(information criteria)\n\t* AIC (Akaike information criterion)\n\t* DIC (Deviance information criterion)\n\t* WAIC (Widely applicable information criterion)\n\t* LOO-CV (Leave-one-out cross-validation)\n\n* 模型平均(model averaging)\n","metadata":{"id":"0B7776895FEF442D85C9D8EB862C5DDE","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"## Example of Model Comparison","metadata":{"id":"B9AE8C90120544888A3C667F5C0668F9","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"![Image Name](https://cdn.kesci.com/upload/image/rkz1ehen1l.png?imageView2/0/w/960/h/960)","metadata":{"id":"856B6574DF224D819DC6DC8437AD8CF4","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"### (1) 提出研究问题\n\nStroop任务是心理学中最常见的认知任务之一，可用于反应抑制能力。\n\n![](https://www.researchgate.net/profile/Ata_Akin/publication/281167153/figure/download/fig1/AS:391418049777669@1470332743733/Three-different-stimulus-conditions-in-the-Stroop-task-neutral-congruent-and.png?_sg=ibeklp8QZ2sbyR29ZZxbOgfS--_RjcKP_uVY36qBahzEJlnMLYPxQyzgYT2Au85eDBClhLqol0A)\n\n个体对于不一致(incongruent)刺激条件下的表现比一致(incongruent)条件下更差，表现为：（1）反应时间往往更长；（2）正确率往往更低。\n\n如果估计不同条件下的个体在反应时间上表现的差异？\n\n图片来源见[这里](https://www.researchgate.net/publication/281167153_Similarity_analysis_of_functional_connectivity_with_functional_near-infrared_spectroscopy/figures?lo=1&utm_source=bing&utm_medium=organic)。","metadata":{"id":"E778994DB1134F4D89776C25EEF221B1","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"### (2) 查看原始数据","metadata":{"id":"6B997D53C7444954BE935657972A85BC","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport arviz as az\nimport pymc3 as pm\nimport matplotlib.pyplot as plt\n\n# theano.tensor是一个对向量进行操作的模块\nimport theano.tensor as tt\n\n# random是一个用于进行随机数操作的包\nimport random","metadata":{"id":"69BBC3E98BD8448F88915E0663F47E44","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"WARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\n"}],"execution_count":1},{"cell_type":"code","source":"raw_data = pd.read_csv('/home/mw/input/data4676/stroop.csv') # 载入数据\nraw_data.head(5) # 查看前5行数据","metadata":{"id":"86E6C9B8E67849C38131C4DFE55081C5","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"        Unnamed: 0             battery_name    condition  correct  \\\n0  stroop_s000_004  Self Regulation Battery  incongruent      0.0   \n1  stroop_s000_006  Self Regulation Battery    congruent      1.0   \n2  stroop_s000_008  Self Regulation Battery    congruent      1.0   \n3  stroop_s000_010  Self Regulation Battery    congruent      1.0   \n4  stroop_s000_012  Self Regulation Battery    congruent      1.0   \n\n   correct_response exp_stage experiment_exp_id           finishtime  \\\n0              66.0  practice            stroop  2016-07-22 04:44:22   \n1              82.0  practice            stroop  2016-07-22 04:44:22   \n2              66.0  practice            stroop  2016-07-22 04:44:22   \n3              71.0  practice            stroop  2016-07-22 04:44:22   \n4              71.0  practice            stroop  2016-07-22 04:44:22   \n\n   focus_shifts  full_screen  key_press possible_responses    rt stim_color  \\\n0             0         True       71.0       [66, 71, 82]  1141       blue   \n1             0         True       82.0       [66, 71, 82]  1048        red   \n2             0         True       66.0       [66, 71, 82]   855       blue   \n3             0         True       71.0       [66, 71, 82]   421      green   \n4             0         True       71.0       [66, 71, 82]   694      green   \n\n  stim_word                                           stimulus  time_elapsed  \\\n0       red  <div class = centerbox><div class = stroop-sti...         24318   \n1       red  <div class = centerbox><div class = stroop-sti...         27574   \n2      blue  <div class = centerbox><div class = stroop-sti...         30828   \n3     green  <div class = centerbox><div class = stroop-sti...         34083   \n4     green  <div class = centerbox><div class = stroop-sti...         37340   \n\n  trial_id           trial_type worker_id  \n0     stim  poldrack-categorize      s001  \n1     stim  poldrack-categorize      s001  \n2     stim  poldrack-categorize      s001  \n3     stim  poldrack-categorize      s001  \n4     stim  poldrack-categorize      s001  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Unnamed: 0</th>\n      <th>battery_name</th>\n      <th>condition</th>\n      <th>correct</th>\n      <th>correct_response</th>\n      <th>exp_stage</th>\n      <th>experiment_exp_id</th>\n      <th>finishtime</th>\n      <th>focus_shifts</th>\n      <th>full_screen</th>\n      <th>key_press</th>\n      <th>possible_responses</th>\n      <th>rt</th>\n      <th>stim_color</th>\n      <th>stim_word</th>\n      <th>stimulus</th>\n      <th>time_elapsed</th>\n      <th>trial_id</th>\n      <th>trial_type</th>\n      <th>worker_id</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>stroop_s000_004</td>\n      <td>Self Regulation Battery</td>\n      <td>incongruent</td>\n      <td>0.0</td>\n      <td>66.0</td>\n      <td>practice</td>\n      <td>stroop</td>\n      <td>2016-07-22 04:44:22</td>\n      <td>0</td>\n      <td>True</td>\n      <td>71.0</td>\n      <td>[66, 71, 82]</td>\n      <td>1141</td>\n      <td>blue</td>\n      <td>red</td>\n      <td>&lt;div class = centerbox&gt;&lt;div class = stroop-sti...</td>\n      <td>24318</td>\n      <td>stim</td>\n      <td>poldrack-categorize</td>\n      <td>s001</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>stroop_s000_006</td>\n      <td>Self Regulation Battery</td>\n      <td>congruent</td>\n      <td>1.0</td>\n      <td>82.0</td>\n      <td>practice</td>\n      <td>stroop</td>\n      <td>2016-07-22 04:44:22</td>\n      <td>0</td>\n      <td>True</td>\n      <td>82.0</td>\n      <td>[66, 71, 82]</td>\n      <td>1048</td>\n      <td>red</td>\n      <td>red</td>\n      <td>&lt;div class = centerbox&gt;&lt;div class = stroop-sti...</td>\n      <td>27574</td>\n      <td>stim</td>\n      <td>poldrack-categorize</td>\n      <td>s001</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>stroop_s000_008</td>\n      <td>Self Regulation Battery</td>\n      <td>congruent</td>\n      <td>1.0</td>\n      <td>66.0</td>\n      <td>practice</td>\n      <td>stroop</td>\n      <td>2016-07-22 04:44:22</td>\n      <td>0</td>\n      <td>True</td>\n      <td>66.0</td>\n      <td>[66, 71, 82]</td>\n      <td>855</td>\n      <td>blue</td>\n      <td>blue</td>\n      <td>&lt;div class = centerbox&gt;&lt;div class = stroop-sti...</td>\n      <td>30828</td>\n      <td>stim</td>\n      <td>poldrack-categorize</td>\n      <td>s001</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>stroop_s000_010</td>\n      <td>Self Regulation Battery</td>\n      <td>congruent</td>\n      <td>1.0</td>\n      <td>71.0</td>\n      <td>practice</td>\n      <td>stroop</td>\n      <td>2016-07-22 04:44:22</td>\n      <td>0</td>\n      <td>True</td>\n      <td>71.0</td>\n      <td>[66, 71, 82]</td>\n      <td>421</td>\n      <td>green</td>\n      <td>green</td>\n      <td>&lt;div class = centerbox&gt;&lt;div class = stroop-sti...</td>\n      <td>34083</td>\n      <td>stim</td>\n      <td>poldrack-categorize</td>\n      <td>s001</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>stroop_s000_012</td>\n      <td>Self Regulation Battery</td>\n      <td>congruent</td>\n      <td>1.0</td>\n      <td>71.0</td>\n      <td>practice</td>\n      <td>stroop</td>\n      <td>2016-07-22 04:44:22</td>\n      <td>0</td>\n      <td>True</td>\n      <td>71.0</td>\n      <td>[66, 71, 82]</td>\n      <td>694</td>\n      <td>green</td>\n      <td>green</td>\n      <td>&lt;div class = centerbox&gt;&lt;div class = stroop-sti...</td>\n      <td>37340</td>\n      <td>stim</td>\n      <td>poldrack-categorize</td>\n      <td>s001</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":2}],"execution_count":2},{"cell_type":"code","source":"random.seed(123) # 设置随机数种子，使每次随机数取样的结果一致\nsamps = random.sample(list(raw_data['worker_id'].unique()), 1) #随机选取1个被试\ndata = raw_data[(raw_data.worker_id.isin(samps)) & (raw_data.exp_stage == \"test\") & (raw_data.rt > 0)] # 选取数据中随机选取到的10名被试在测试阶段的反应时大于0的数据\ndata = data[[\"worker_id\",\"correct\",\"condition\", \"rt\"]] # 选取数据中的判断正确率，刺激条件，和被试编号\ndata.reset_index(inplace=True, drop=True) # 重置每个试次的编号\ndata.head(5) # 查看数据的前5行","metadata":{"id":"B3299CD3E52C45F8814772C40242490C","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"  worker_id  correct    condition   rt\n0      s058      1.0  incongruent  480\n1      s058      1.0    congruent  576\n2      s058      1.0  incongruent  757\n3      s058      1.0    congruent  760\n4      s058      1.0    congruent  640","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>worker_id</th>\n      <th>correct</th>\n      <th>condition</th>\n      <th>rt</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>s058</td>\n      <td>1.0</td>\n      <td>incongruent</td>\n      <td>480</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>s058</td>\n      <td>1.0</td>\n      <td>congruent</td>\n      <td>576</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>s058</td>\n      <td>1.0</td>\n      <td>incongruent</td>\n      <td>757</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>s058</td>\n      <td>1.0</td>\n      <td>congruent</td>\n      <td>760</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>s058</td>\n      <td>1.0</td>\n      <td>congruent</td>\n      <td>640</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":3}],"execution_count":3},{"cell_type":"code","source":"def filter_ab(df):\n    '''\n    判断数据是否在平均值三个标准差以内\n    '''\n    mp1 = df.mean() - 3 * df.std()\n    mp2 = df.mean() + 3 * df.std()\n    norm = df.where((df > mp1) & (df < mp2)) # 判断数据是否在三个标准差以内，不在则赋值为NA\n    return norm\ndata.rt = data.groupby('worker_id').rt.apply(filter_ab) # 将数据按照被试进行分组，分别检查每个被试的反应时是否在三个标准差以内\ndata = data.dropna(axis=0, how='any') # 将含NAN的行进行删除\ndata.reset_index(inplace=True, drop=True) # 重置每个试次的编号\n","metadata":{"id":"5B3B91E4343E4702B21EF448F49AF6B5","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":4},{"cell_type":"code","source":"# 按照条件进行分组，绘制反应时间的概率密度图\ndata.groupby('condition').rt.describe()","metadata":{"id":"C786BDC467894CD3915387A2C3BA2636","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"             count        mean         std    min     25%    50%     75%  \\\ncondition                                                                  \ncongruent     48.0  602.145833  151.296664  367.0  509.00  552.0  646.25   \nincongruent   48.0  726.729167  151.551794  480.0  608.25  718.0  847.00   \n\n                max  \ncondition            \ncongruent    1144.0  \nincongruent  1066.0  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>count</th>\n      <th>mean</th>\n      <th>std</th>\n      <th>min</th>\n      <th>25%</th>\n      <th>50%</th>\n      <th>75%</th>\n      <th>max</th>\n    </tr>\n    <tr>\n      <th>condition</th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>congruent</th>\n      <td>48.0</td>\n      <td>602.145833</td>\n      <td>151.296664</td>\n      <td>367.0</td>\n      <td>509.00</td>\n      <td>552.0</td>\n      <td>646.25</td>\n      <td>1144.0</td>\n    </tr>\n    <tr>\n      <th>incongruent</th>\n      <td>48.0</td>\n      <td>726.729167</td>\n      <td>151.551794</td>\n      <td>480.0</td>\n      <td>608.25</td>\n      <td>718.0</td>\n      <td>847.00</td>\n      <td>1066.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":5}],"execution_count":5},{"cell_type":"code","source":"# 按照条件进行分组，绘制反应时间的概率密度图\ndata.groupby(['condition']).rt.plot.density() ","metadata":{"id":"6718558B360D444CB62DAB747A2EC744","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"condition\ncongruent      AxesSubplot(0.125,0.125;0.775x0.755)\nincongruent    AxesSubplot(0.125,0.125;0.775x0.755)\nName: rt, dtype: object"},"metadata":{},"execution_count":6},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/6718558B360D444CB62DAB747A2EC744/rloww16zhc.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":6},{"cell_type":"code","source":"# 按照条件进行分组，绘制反应时间的直方图\ndata.groupby('condition').rt.plot.hist(bins=30, alpha=0.5)","metadata":{"id":"643ADE4018844623A4203B21892726EC","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"condition\ncongruent      AxesSubplot(0.125,0.125;0.775x0.755)\nincongruent    AxesSubplot(0.125,0.125;0.775x0.755)\nName: rt, dtype: object"},"metadata":{},"execution_count":7},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/643ADE4018844623A4203B21892726EC/rloww1fhk7.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":7},{"cell_type":"code","source":"# 将‘condition’条件进行0， 1编码，congruent转化为1，incongruent转化为0，以便在pymc3中进行计算\ndata.condition = data.condition.map({'congruent':1, 'incongruent':0})","metadata":{"id":"DAC44E2648A74AEAB77118CC7FC02AFB","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":8},{"cell_type":"markdown","source":"从观察数据中获得的信息：\n* 反应时间数据是“连续”的数据；\n* 反应时间的分布看起来不太像是正常分布；\n* 两种条件下的差异比较微妙。\n* 不同被试的反应时分布存在差异","metadata":{"id":"341F7FA5C9364EBCAB99062EF7CB68E6","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"可能的模型：\n* 常规做法，以高低冲突作为自变量，反应时间作为因变量，采用线性模型（即正态分布作为likelihood），检验高低冲突对反应时总体的均值的影响；\n* 广义线性模型，尝试非正态分布作为likelihood：log-normal，gamma, ex-Gausian.\n\n通过贝叶斯分析来比较来推断高低冲突对反应时间的影响，其中，需要进行模型比较来选择最优的模型。","metadata":{"id":"3BD8645007B04324B91D2206921353E7","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"#### Model 1: Normal\n\n首先，我们假定反应时是呈正态分布的，正态分布的均值是由高一致性条件和低一致性条件决定的。\n\n\n![Image Name](https://docs.pymc.io/en/v3/_images/continuous-18.png)\n\n图片来源：https://docs.pymc.io/en/v3/api/distributions/continuous.html#pymc3.distributions.continuous.Normal","metadata":{"id":"EC48D74A214F4690BEF49CA228092816","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 模型建构与先验选择\n\n写作线性模型的概率形式：\n\n$y_i = \\alpha + \\beta*x_i + e$\n\n$y_i \\sim N(\\mu, \\sigma)$\n\n$\\mu = \\alpha + \\beta * cond_{i}$\n\n$\\alpha \\sim N(700, 200)$\n\n$\\beta \\sim N(100, 200)$\n\n$\\sigma \\sim HN(100)$","metadata":{"id":"BB494F96A3944B5789E7B951EEB1C5BE","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with pm.Model() as NormalModel:\n    # 先验分布: alpha, beta, sigma这三个参数是随机变量\n    sigma = pm.HalfNormal('sigma', sd=200)\n    alpha = pm.Normal('alpha', mu=700, sd=200)\n    beta = pm.Normal('beta', mu=100, sd=200)\n    # 自变量condition是之前已经载入的数据\n    x = pm.Data(\"x\", data['condition'])\n    # 正态分布均值是确定性随机变量，这个变量的值完全由右端值确定\n    mu = pm.Deterministic(\"mu\", alpha + beta*x) \n    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    # 假定因变量服从正态分布\n    y_obs = pm.Normal('y_obs', mu=mu, sd=sigma, observed=data['rt'] )","metadata":{"id":"856AB9E42016419E8D5139F6477C42E5","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":9},{"cell_type":"code","source":"# 展示模型结构\npm.model_to_graphviz(NormalModel)","metadata":{"id":"2D0703F0B4AB4F4C846D2E27355DE429","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7f30a744d150>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/2D0703F0B4AB4F4C846D2E27355DE429/rloww4yzgs.svg\">"},"metadata":{},"execution_count":10}],"execution_count":10},{"cell_type":"markdown","source":"##### 先验检验(prior predictive check)\n\n先验的设定是否合理？\n\n可以通过先验预测检验（ Prior Predictive Distribution ）来进行初步的判断。\n\n先验预测检验：利用模型和先验生成假数据并利用这些生成的假数据来评估先验是否合理的过程。","metadata":{"id":"8E07402F752B46DF82BD1932976AE40D","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with NormalModel:\n    # 先验预测检查\n    prior_checks = pm.sample_prior_predictive(samples=50)","metadata":{"id":"994B75EFB1584A088C92A7383183CA32","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":11},{"cell_type":"code","source":"x = np.random.randint(2, size = 50) #生成50个假数据，取值为[0,1]\n\nfor a, b in zip(prior_checks[\"alpha\"], prior_checks[\"beta\"]):\n    y = a + b * x \n    plt.plot(x, y)","metadata":{"id":"682E96A1EA25453F905C820E8EB2461E","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/682E96A1EA25453F905C820E8EB2461E/rloww56vt1.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":12},{"cell_type":"markdown","source":"通过先验预测检查，我们发现一致和不一致条件的因变量的取值分别在200-1200, 0-1400之间，还算合理。","metadata":{"id":"D1B83C1368B346F69D89535BB11F7CED","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 计算后验","metadata":{"id":"67F99683715244188C6DB265A604B212","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with NormalModel:\n    # 不直接使用InferenceData的理由是为了便于之后的模型比较\n    trace1_for_comp = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2)    \n    # 将pymc的采样对象转化为inferencedata\n    trace1=az.from_pymc3(trace1_for_comp)","metadata":{"id":"FF56A7C41EC446D697795DEFABE8579B","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:3: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.\n  This is separate from the ipykernel package so we can avoid doing imports until\nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [beta, alpha, sigma]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 6 seconds.\n"}],"execution_count":13},{"cell_type":"markdown","source":"##### MCMC诊断\n\n计算后验分布之后，我们需要检查MCMC的采样情况是否稳定，我们一般采取目视检查法和统计数据$\\hat{R}$。","metadata":{"id":"28EEF0F1F7C54BEF86158E7F3E1ECA96","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"# 绘制各参数的采样情况\naz.plot_trace(trace1,var_names=['alpha','beta','sigma'])","metadata":{"id":"3CCB9BE932F44FF48E549072BE20CC10","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"array([[<AxesSubplot:title={'center':'alpha'}>,\n        <AxesSubplot:title={'center':'alpha'}>],\n       [<AxesSubplot:title={'center':'beta'}>,\n        <AxesSubplot:title={'center':'beta'}>],\n       [<AxesSubplot:title={'center':'sigma'}>,\n        <AxesSubplot:title={'center':'sigma'}>]], dtype=object)"},"metadata":{},"execution_count":14},{"output_type":"display_data","data":{"text/plain":"<Figure size 864x432 with 6 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/3CCB9BE932F44FF48E549072BE20CC10/rlowwg4qks.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":14},{"cell_type":"code","source":"az.summary(trace1, var_names=[\"alpha\", \"beta\", \"sigma\"])","metadata":{"id":"E4F937E679B84AC6830D9DFD526BF3A7","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"          mean      sd   hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\nalpha  723.592  22.477  682.228  767.571      0.535    0.379    1766.0   \nbeta  -118.837  31.347 -176.663  -60.538      0.775    0.548    1640.0   \nsigma  152.954  11.371  132.294  173.876      0.236    0.167    2303.0   \n\n       ess_tail  r_hat  \nalpha    2211.0    1.0  \nbeta     2539.0    1.0  \nsigma    2439.0    1.0  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>alpha</th>\n      <td>723.592</td>\n      <td>22.477</td>\n      <td>682.228</td>\n      <td>767.571</td>\n      <td>0.535</td>\n      <td>0.379</td>\n      <td>1766.0</td>\n      <td>2211.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>beta</th>\n      <td>-118.837</td>\n      <td>31.347</td>\n      <td>-176.663</td>\n      <td>-60.538</td>\n      <td>0.775</td>\n      <td>0.548</td>\n      <td>1640.0</td>\n      <td>2539.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>sigma</th>\n      <td>152.954</td>\n      <td>11.371</td>\n      <td>132.294</td>\n      <td>173.876</td>\n      <td>0.236</td>\n      <td>0.167</td>\n      <td>2303.0</td>\n      <td>2439.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":15}],"execution_count":15},{"cell_type":"markdown","source":"根据采样得到的trace比较稳定，并且r_hat取值接近1说明MCMC采样结果稳定。","metadata":{"id":"431661E653E84EB08E4469FAA108F87B","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 模型评估(PPC)\n\n在该模型的采样情况达到稳定之后，我们需要判断该模型能够预测真实数据，这就需要使用后验预测分布进行检查。","metadata":{"id":"E04D6EDE72654968AEB1ABA55FCCF83C","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with NormalModel:\n     # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n    ppc_y = pm.sample_posterior_predictive(trace1.posterior) \n#将ppc_y转化为InferenceData对象合并到trace中\naz.concat(trace1, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)","metadata":{"id":"C8979DB6A7BD4C94BF5A1482111F1765","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"update_display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n    <div>\n      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      100.00% [4000/4000 00:05&lt;00:00]\n    </div>\n    "},"metadata":{}},{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n"}],"execution_count":16},{"cell_type":"code","source":"# 绘制后验预测分布\naz.plot_ppc(trace1)","metadata":{"id":"233B1628048942568450321282BF4F53","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<AxesSubplot:xlabel='y_obs'>"},"metadata":{},"execution_count":17},{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/IPython/core/pylabtools.py:151: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n  fig.canvas.print_figure(bytes_io, **kw)\n"},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/233B1628048942568450321282BF4F53/rlowxas7rc.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":17},{"cell_type":"markdown","source":"#### Model2: Log-normal\n\n然后，我们假定反应时是呈log正态分布的，LogNormal分布的均值是由高一致性条件和低一致性条件决定的。\n\n![Image Name](https://docs.pymc.io/en/v3/_images/continuous-16.png)\n\n\n图片来源：https://docs.pymc.io/en/v3/api/distributions/continuous.html#pymc3.distributions.continuous.LogitNormal","metadata":{"id":"B0570E5F803540CEB75D0508044F2F86","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 模型建构与先验选择\n\n写作线性模型的概率形式：\n\n$y_i \\sim lognormal(\\mu, \\sigma)$   # 将normal 换成lognormal\n\n$\\mu = \\alpha + \\beta * cond_{i}$\n\n$\\alpha \\sim N(700, 200)$\n\n$\\beta \\sim N(100, 200)$\n\n$\\sigma \\sim HN(100)$","metadata":{"id":"41CB695EE37F4EC7874D66F3BD9AF98B","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with pm.Model() as LogNormal:\n    # 先验分布: alpha, beta, sigma这三个参数是随机变量\n    sigma = pm.HalfNormal('sigma', sd=100)\n    alpha = pm.Normal('alpha', mu=700, sd=200)\n    beta = pm.Normal('beta', mu=100, sd=200)\n    # 自变量conf是之前已经载入的数据\n    x = pm.Data(\"x\", data['condition'])\n    # 正态分布均值是确定性随机变量，这个变量的值完全由右端值确定\n    mu = pm.Deterministic(\"mu\", alpha + beta*x) \n    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    # 假定因变量服从log-normal分布\n    y_obs = pm.Lognormal('y_obs',mu=mu,sd=sigma,observed=data['rt'] )","metadata":{"id":"4541DCBB9FD740E68288C6D3EE7597ED","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":18},{"cell_type":"code","source":"# 展示模型结构\npm.model_to_graphviz(LogNormal)","metadata":{"id":"924354FC9A7F44AA976E0960DFD9BAF4","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7f30a598e310>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/924354FC9A7F44AA976E0960DFD9BAF4/rlowxa3jo5.svg\">"},"metadata":{},"execution_count":19}],"execution_count":19},{"cell_type":"markdown","source":"##### 先验检验(prior predictive check)\n\n先验的设定是否合理？\n\n可以通过先验预测检验（ Prior Predictive Distribution ）来进行初步的判断。\n\n先验预测检验：利用模型和先验生成假数据并利用这些生成的假数据来评估先验是否合理的过程。","metadata":{"id":"C3406EF75EEC4CAEAFB568C559F63F0E","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with LogNormal:\n    # 先验预测检查\n    prior_checks = pm.sample_prior_predictive(samples=50)","metadata":{"id":"92FC4C97245F47459FD911BD9A9AC5AB","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/pymc3/distributions/continuous.py:1864: RuntimeWarning: overflow encountered in exp\n  return np.exp(mu + (tau ** -0.5) * samples)\n"}],"execution_count":20},{"cell_type":"code","source":"x = np.random.randint(2, size = 50) #生成50个假数据，取值为[0,1]\n\nfor a, b in zip(prior_checks[\"alpha\"], prior_checks[\"beta\"]):\n    y = a + b * x \n    plt.plot(x, y)","metadata":{"id":"80E911BC84CF4E0EA69CA5DFE8685EBE","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/80E911BC84CF4E0EA69CA5DFE8685EBE/rlowxbojrk.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":21},{"cell_type":"markdown","source":"通过先验预测检查，我们发现一致和不一致条件的因变量的取值分别在400-1000，200-1400之间，比较合理。","metadata":{"id":"9F697FB7B6634DBD8700979C27F45B2C","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 计算后验","metadata":{"id":"02E83865767A473E8D68ED96F61772D3","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with LogNormal:\n    # 不直接使用InferenceData的理由是为了便于之后的模型比较\n    trace2_for_comp = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2)    \n    # 将pymc的采样对象转化为inferencedata\n    trace2=az.from_pymc3(trace2_for_comp)","metadata":{"id":"E29CE71649724D05A112D15701A93128","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:3: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.\n  This is separate from the ipykernel package so we can avoid doing imports until\nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [beta, alpha, sigma]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 9 seconds.\n"}],"execution_count":22},{"cell_type":"markdown","source":"##### MCMC诊断\n\n计算后验分布之后，我们需要检查MCMC的采样情况是否稳定，我们一般采取目视检查法和统计数据r_hat","metadata":{"id":"8F3192A3561E469496D146C59294B330","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"# 绘制各参数的采样情况\naz.plot_trace(trace1,var_names=['alpha','beta','sigma'])","metadata":{"id":"55028FE26D58425BA7CC139D2326648E","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"array([[<AxesSubplot:title={'center':'alpha'}>,\n        <AxesSubplot:title={'center':'alpha'}>],\n       [<AxesSubplot:title={'center':'beta'}>,\n        <AxesSubplot:title={'center':'beta'}>],\n       [<AxesSubplot:title={'center':'sigma'}>,\n        <AxesSubplot:title={'center':'sigma'}>]], dtype=object)"},"metadata":{},"execution_count":23},{"output_type":"display_data","data":{"text/plain":"<Figure size 864x432 with 6 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/55028FE26D58425BA7CC139D2326648E/rlowxoig1.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":23},{"cell_type":"code","source":"az.summary(trace2, var_names=[\"alpha\", \"beta\", \"sigma\"])","metadata":{"id":"4C823DE16583474A9B4571BABDFD33D4","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"        mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\nalpha  6.567  0.032   6.505    6.625      0.001    0.001    1765.0    2183.0   \nbeta  -0.193  0.046  -0.277   -0.106      0.001    0.001    1901.0    2097.0   \nsigma  0.219  0.017   0.189    0.251      0.000    0.000    2214.0    2490.0   \n\n       r_hat  \nalpha    1.0  \nbeta     1.0  \nsigma    1.0  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>alpha</th>\n      <td>6.567</td>\n      <td>0.032</td>\n      <td>6.505</td>\n      <td>6.625</td>\n      <td>0.001</td>\n      <td>0.001</td>\n      <td>1765.0</td>\n      <td>2183.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>beta</th>\n      <td>-0.193</td>\n      <td>0.046</td>\n      <td>-0.277</td>\n      <td>-0.106</td>\n      <td>0.001</td>\n      <td>0.001</td>\n      <td>1901.0</td>\n      <td>2097.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>sigma</th>\n      <td>0.219</td>\n      <td>0.017</td>\n      <td>0.189</td>\n      <td>0.251</td>\n      <td>0.000</td>\n      <td>0.000</td>\n      <td>2214.0</td>\n      <td>2490.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":24}],"execution_count":24},{"cell_type":"markdown","source":"根据采样得到的trace比较稳定，并且$\\hat{R}$取值接近1说明MCMC采样结果稳定。","metadata":{"id":"73F15A16C2014D12BDCB1C10D0228234","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 模型评估(PPC)\n\n在该模型的采样情况达到稳定之后，我们需要判断该模型能够预测真实数据，这就需要使用后验预测分布进行检查。","metadata":{"id":"8C22E820544D45778269188198B8617A","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with LogNormal:\n     # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n    ppc_y = pm.sample_posterior_predictive(trace2.posterior) \n#将ppc_y转化为InferenceData对象合并到trace中\naz.concat(trace2, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)","metadata":{"id":"A635CE6012D04F0DAFC210589FDA4700","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"update_display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n    <div>\n      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      100.00% [4000/4000 00:05&lt;00:00]\n    </div>\n    "},"metadata":{}},{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n"}],"execution_count":25},{"cell_type":"code","source":"# 绘制后验预测分布\naz.plot_ppc(trace2)","metadata":{"id":"35A921DA2384401C894FA8D2850FF304","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<AxesSubplot:xlabel='y_obs'>"},"metadata":{},"execution_count":26},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/35A921DA2384401C894FA8D2850FF304/rlowyf18xh.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":26},{"cell_type":"markdown","source":"#### Model3: Gamma\n\n我们假定反应时是呈Gamma分布。\n\nGamma分布有两个参数$\\alpha$, $\\beta$，类似于正态分布中的的$\\mu$, $\\sigma$。\n\n检验$\\alpha$是否受到冲突水平的影响。\n\n\n![Image Name](https://docs.pymc.io/en/v3/_images/continuous-6.png)\n\n图片来源：https://docs.pymc.io/en/v3/api/distributions/continuous.html#pymc3.distributions.continuous.Gamma","metadata":{"id":"42FDB30F030E4064BABF0F12DC188342","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 模型建构与先验选择\n\n写作线性模型的概率形式：\n\n$y_i \\sim gamma(\\mu, \\sigma)$   # 将normal 换成gamma\n\n$\\mu = \\alpha + \\beta * cond_{i}$\n\n$\\alpha \\sim N(700, 200)$\n\n$\\beta \\sim N(100, 200)$\n\n$\\sigma \\sim HN(100)$\n","metadata":{"id":"2C107B757D634A1FA41655A25AABDC7A","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with pm.Model() as Gamma:\n    # 先验分布: alpha, beta, sigma这三个参数是随机变量\n    sigma = pm.HalfNormal('sigma', sd=100)\n    alpha = pm.Normal('alpha', mu=700, sd=200)\n    beta = pm.Normal('beta', mu=100, sd=200)\n    # 自变量conf是之前已经载入的数据\n    x = pm.Data(\"x\", data['condition'])\n    # 正态分布均值是确定性随机变量，这个变量的值完全由右端值确定\n    mu = pm.Deterministic(\"mu\", alpha + beta*x) \n    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    # 假定因变量服从log-normal分布\n    y_obs = pm.Gamma('y_obs',mu=mu,sd=sigma,observed=data['rt'] )","metadata":{"id":"03B7E63F5AEC490AAB4FF6EAB06AAA63","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":27},{"cell_type":"code","source":"# 展示模型结构\npm.model_to_graphviz(Gamma)","metadata":{"id":"9F53050CB40C429F827242C4251A2FB3","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7f2f6a37cad0>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/9F53050CB40C429F827242C4251A2FB3/rlowyfo5wt.svg\">"},"metadata":{},"execution_count":28}],"execution_count":28},{"cell_type":"markdown","source":"##### 先验检验(prior predictive check)\n\n先验的设定是否合理？\n\n可以通过先验预测检验（ Prior Predictive Distribution ）来进行初步的判断。\n\n先验预测检验：利用模型和先验生成假数据并利用这些生成的假数据来评估先验是否合理的过程。","metadata":{"id":"7358052588AE4FBA936C3484859E40F4","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with Gamma:\n    # 先验预测检查\n    prior_checks = pm.sample_prior_predictive(samples=50)","metadata":{"id":"8148AC1D8A054AC081B5114150184B1B","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":29},{"cell_type":"code","source":"x = np.random.randint(2, size = 50) #生成50个假数据，取值为[0,1]\n\nfor a, b in zip(prior_checks[\"alpha\"], prior_checks[\"beta\"]):\n    y = a + b * x \n    plt.plot(x, y)","metadata":{"id":"44DD5C427DBB40B09D9399631E2BBF55","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/44DD5C427DBB40B09D9399631E2BBF55/rlowygnzg9.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":30},{"cell_type":"markdown","source":"通过先验预测检查，我们发现一致和不一致条件的因变量的取值分别在< 200 ms ~ 1200 ms，400 - 1300之间，比较合理","metadata":{"id":"3E220DAE095E4EE189B8B62511901DE6","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 计算后验","metadata":{"id":"298F9EE466F441A497A64C6D66D64C54","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with Gamma:\n    # 不直接使用InferenceData的理由是为了便于之后的模型比较\n    trace3_for_comp = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2)    \n    # 将pymc的采样对象转化为inferencedata\n    trace3 = az.from_pymc3(trace3_for_comp)","metadata":{"id":"5F6956318E984BE682D83E0EB4EE6725","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:3: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.\n  This is separate from the ipykernel package so we can avoid doing imports until\nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [beta, alpha, sigma]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 6 seconds.\n"}],"execution_count":31},{"cell_type":"markdown","source":"##### MCMC诊断\n\n计算后验分布之后，我们需要检查MCMC的采样情况是否稳定，我们一般采取目视检查法和统计数据r_hat","metadata":{"id":"55AB55A4F03F40978D4429C49DCC5824","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"# 绘制各参数的采样情况\naz.plot_trace(trace3,var_names=['alpha','beta','sigma'])","metadata":{"id":"CD914DAA85E74072861215F3D65C46D4","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"array([[<AxesSubplot:title={'center':'alpha'}>,\n        <AxesSubplot:title={'center':'alpha'}>],\n       [<AxesSubplot:title={'center':'beta'}>,\n        <AxesSubplot:title={'center':'beta'}>],\n       [<AxesSubplot:title={'center':'sigma'}>,\n        <AxesSubplot:title={'center':'sigma'}>]], dtype=object)"},"metadata":{},"execution_count":32},{"output_type":"display_data","data":{"text/plain":"<Figure size 864x432 with 6 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/CD914DAA85E74072861215F3D65C46D4/rlowytum60.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":32},{"cell_type":"code","source":"az.summary(trace3, var_names=[\"alpha\", \"beta\", \"sigma\"])","metadata":{"id":"FC825D01B9F54B828E8C43DA89D8A3EB","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"          mean      sd   hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\nalpha  720.921  20.915  680.127  758.871      0.494    0.349    1789.0   \nbeta  -112.826  28.578 -168.594  -62.618      0.635    0.449    2016.0   \nsigma  145.374  11.008  125.011  165.154      0.223    0.158    2483.0   \n\n       ess_tail  r_hat  \nalpha    1792.0    1.0  \nbeta     1988.0    1.0  \nsigma    2552.0    1.0  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>alpha</th>\n      <td>720.921</td>\n      <td>20.915</td>\n      <td>680.127</td>\n      <td>758.871</td>\n      <td>0.494</td>\n      <td>0.349</td>\n      <td>1789.0</td>\n      <td>1792.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>beta</th>\n      <td>-112.826</td>\n      <td>28.578</td>\n      <td>-168.594</td>\n      <td>-62.618</td>\n      <td>0.635</td>\n      <td>0.449</td>\n      <td>2016.0</td>\n      <td>1988.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>sigma</th>\n      <td>145.374</td>\n      <td>11.008</td>\n      <td>125.011</td>\n      <td>165.154</td>\n      <td>0.223</td>\n      <td>0.158</td>\n      <td>2483.0</td>\n      <td>2552.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":33}],"execution_count":33},{"cell_type":"markdown","source":"根据采样得到的trace比较稳定，并且r_hat取值接近1说明MCMC采样结果稳定。","metadata":{"id":"1688022D64A44BC6A2C3423EB17AEDAD","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 模型评估(PPC)\n\n在该模型的采样情况达到稳定之后，我们需要判断该模型能够预测真实数据，这就需要使用后验预测分布进行检查。","metadata":{"id":"5B2013EBDBF6406E873CC2F1D80527FF","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with Gamma:\n     # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n    ppc_y = pm.sample_posterior_predictive(trace3.posterior) \n#将ppc_y转化为InferenceData对象合并到trace中\naz.concat(trace3, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)","metadata":{"id":"0A2FBB20510141C591B6D71A2E6AA92D","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"update_display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n    <div>\n      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      100.00% [4000/4000 00:06&lt;00:00]\n    </div>\n    "},"metadata":{}},{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n"}],"execution_count":34},{"cell_type":"code","source":"# 绘制后验预测分布\naz.plot_ppc(trace3)","metadata":{"id":"CC7AA857F000425483AE183C7C5514B2","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<AxesSubplot:xlabel='y_obs'>"},"metadata":{},"execution_count":35},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/CC7AA857F000425483AE183C7C5514B2/rlowzhz4sa.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":35},{"cell_type":"markdown","source":"#### Model4: ex-Gaussian\n\n我们假定反应时是ex-Gaussian 模型。\n\n该模型有三个参数：$\\mu$, $\\sigma$, $\\tau$\n\n![Image Name](https://docs.pymc.io/en/latest/_images/pymc-ExGaussian-1.png)\n\n图片来源：https://docs.pymc.io/en/latest/api/distributions/generated/pymc.ExGaussian.html","metadata":{"id":"12A9F25D44D84305824A8040EE2D2251","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 模型建构与先验选择\n\n写作线性模型的概率形式：\n\n$y_i \\sim exgaussian(\\mu, \\sigma, \\tau)$   # 将normal 换成 exgaussian\n\n$\\mu = \\alpha + \\beta * cond_{i}$\n\n$\\alpha \\sim N(700,200)$\n\n$\\beta \\sim N(100,200)$\n\n$\\sigma \\sim HN(100)$\n\n$\\tau \\sim HN(100)$","metadata":{"id":"914B8918639249C48E2717208E520486","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with pm.Model() as exgaussian:\n    # 先验分布:alpha,beta,sigma,nu这三个参数是随机变量\n    alpha = pm.Normal('alpha', mu = 700, sd=200)\n    beta = pm.Normal('beta',mu=100, sd=200)\n    sigma = pm.HalfNormal('sigma', sd=200) \n    nu = pm.HalfNormal('nu', sd=100)   \n    # 自变量conf是之前已经载入的数据\n    x = pm.Data(\"x\", data['condition'])\n    # 参数mu是确定性随机变量，这个变量的值完全由右端值确定\n    mu = pm.Deterministic(\"mu\",  alpha+ beta*x) \n\n    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    y_obs = pm.ExGaussian('y_obs',mu=mu,sigma=sigma,nu=nu,observed=data['rt'] )","metadata":{"id":"DEA4E7C2DCDC4EEFA77D50FA56B86FEE","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":36},{"cell_type":"code","source":"pm.model_to_graphviz(exgaussian)","metadata":{"id":"3A4F883EBD5D4560B490B19E44B03310","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7f2f66a4b1d0>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/3A4F883EBD5D4560B490B19E44B03310/rlowzhtt54.svg\">"},"metadata":{},"execution_count":37}],"execution_count":37},{"cell_type":"markdown","source":"##### 先验检验(prior predictive check)\n\n先验的设定是否合理？\n\n可以通过先验预测检验（ Prior Predictive Distribution ）来进行初步的判断。\n\n先验预测检验：利用模型和先验生成假数据并利用这些生成的假数据来评估先验是否合理的过程。","metadata":{"id":"10ADE0FF931448B5857AB2B8EE6F7241","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with exgaussian:\n    \n    # 先验预测检查\n    prior_checks = pm.sample_prior_predictive(samples=50)","metadata":{"id":"8D915F3272A548EF87DB4B9AEBC8F9A2","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":38},{"cell_type":"code","source":"x = np.random.randint(2, size = 50) #生成50个假数据，取值为[0,1]\n\nfor a, b in zip(prior_checks[\"alpha\"], prior_checks[\"beta\"]):\n    y = a + b * x \n    plt.plot(x, y)","metadata":{"id":"965E309247724FE3A0A5CEFD2DDD76AF","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/965E309247724FE3A0A5CEFD2DDD76AF/rlowzid9rd.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":39},{"cell_type":"markdown","source":"##### 计算后验","metadata":{"id":"772ACCC8FBE44D1BBD4AFF6AE48D6619","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with exgaussian:\n    trace4_for_comp = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2)    \n    # 将pymc的采样对象转化为inferencedata\n    trace4=az.from_pymc3(trace4_for_comp)","metadata":{"id":"A3DF23DB2FDB46AA82A1823264B537FE","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"WARNING (theano.tensor.opt): Optimization Warning: The Op erfcx does not provide a C implementation. As well as being potentially slow, this also disables loop fusion.\n/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:2: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.\n  \nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [nu, sigma, beta, alpha]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 13 seconds.\n"}],"execution_count":40},{"cell_type":"markdown","source":"##### MCMC诊断\n\n计算后验分布之后，我们需要检查MCMC的采样情况是否稳定，我们一般采取目视检查法和统计数据r_hat","metadata":{"id":"9ABF31DD0FC444D780E0326A2FB390CD","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"az.summary(trace4, var_names=[\"alpha\", \"beta\", \"sigma\",\"nu\"])","metadata":{"id":"BA9205F9CB844F4E9B013469A7CF8E7D","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"          mean      sd   hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\nalpha  572.880  30.742  515.869  627.060      0.922    0.654    1120.0   \nbeta   -98.337  26.982 -148.513  -48.113      0.728    0.515    1373.0   \nsigma   71.067  19.261   37.898  106.843      0.535    0.382    1318.0   \nnu     140.644  25.369   94.483  188.830      0.676    0.478    1425.0   \n\n       ess_tail  r_hat  \nalpha    1389.0    1.0  \nbeta     2001.0    1.0  \nsigma    1830.0    1.0  \nnu       1407.0    1.0  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>alpha</th>\n      <td>572.880</td>\n      <td>30.742</td>\n      <td>515.869</td>\n      <td>627.060</td>\n      <td>0.922</td>\n      <td>0.654</td>\n      <td>1120.0</td>\n      <td>1389.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>beta</th>\n      <td>-98.337</td>\n      <td>26.982</td>\n      <td>-148.513</td>\n      <td>-48.113</td>\n      <td>0.728</td>\n      <td>0.515</td>\n      <td>1373.0</td>\n      <td>2001.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>sigma</th>\n      <td>71.067</td>\n      <td>19.261</td>\n      <td>37.898</td>\n      <td>106.843</td>\n      <td>0.535</td>\n      <td>0.382</td>\n      <td>1318.0</td>\n      <td>1830.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>nu</th>\n      <td>140.644</td>\n      <td>25.369</td>\n      <td>94.483</td>\n      <td>188.830</td>\n      <td>0.676</td>\n      <td>0.478</td>\n      <td>1425.0</td>\n      <td>1407.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":41}],"execution_count":41},{"cell_type":"code","source":"# 绘制各参数的采样情况\naz.plot_trace(trace4,var_names=[\"alpha\", \"beta\", \"sigma\",\"nu\"])","metadata":{"id":"6BCA8623398B483BBCDF67E1ADF28FA1","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"array([[<AxesSubplot:title={'center':'alpha'}>,\n        <AxesSubplot:title={'center':'alpha'}>],\n       [<AxesSubplot:title={'center':'beta'}>,\n        <AxesSubplot:title={'center':'beta'}>],\n       [<AxesSubplot:title={'center':'sigma'}>,\n        <AxesSubplot:title={'center':'sigma'}>],\n       [<AxesSubplot:title={'center':'nu'}>,\n        <AxesSubplot:title={'center':'nu'}>]], dtype=object)"},"metadata":{},"execution_count":42},{"output_type":"display_data","data":{"text/plain":"<Figure size 864x576 with 8 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/6BCA8623398B483BBCDF67E1ADF28FA1/rlox03kexf.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":42},{"cell_type":"markdown","source":"根据采样得到的trace比较稳定，并且r_hat取值接近1说明MCMC采样结果稳定。","metadata":{"id":"2D46F96863EB407683308D6EE69AA612","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"##### 模型评估(PPC)\n\n在该模型的采样情况达到稳定之后，我们需要判断该模型能够预测真实数据，这就需要使用后验预测分布进行检查。","metadata":{"id":"2112A62BF08640F88E6736590DFD9130","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with exgaussian:\n    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n    ppc_y = pm.sample_posterior_predictive(trace4.posterior) \n#将ppc_y转化为InferenceData对象合并到trace中\naz.concat(trace4, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)","metadata":{"id":"9FE44B116D6646848546867295D76230","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"update_display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n    <div>\n      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      100.00% [4000/4000 00:03&lt;00:00]\n    </div>\n    "},"metadata":{}},{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n"}],"execution_count":43},{"cell_type":"code","source":"az.plot_ppc(trace4)","metadata":{"id":"4B4781A01D944D1E889DDE114F650BC9","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<AxesSubplot:xlabel='y_obs'>"},"metadata":{},"execution_count":44},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/4B4781A01D944D1E889DDE114F650BC9/rlox0veq9w.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":44},{"cell_type":"markdown","source":"### 模型比较\n\n#### LOO","metadata":{"id":"7CFD134A469545599499E9FA74DD80F7","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"# 将四个模型的采样结果进行比较\ncompare_dict = {\"normal\": trace1, \"log-normal\": trace2, \"gamma\": trace3,\"exgaussian\":trace4}\n# 选择loo方法进行比较\ncomp = az.compare(compare_dict, ic='loo')\ncomp","metadata":{"id":"CCE5C255DDCE4A5B8F4F5BD97E1ED386","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n  \"The default method used to estimate the weights for each model,\"\n"},{"output_type":"execute_result","data":{"text/plain":"            rank         loo     p_loo     d_loo        weight        se  \\\nlog-normal     0 -612.674051  3.196749  0.000000  6.360410e-01  8.412337   \nexgaussian     1 -613.321656  4.551529  0.647606  3.639590e-01  7.838731   \ngamma          2 -614.527448  2.904233  1.853397  0.000000e+00  7.928021   \nnormal         3 -620.535490  3.511112  7.861439  1.287859e-14  8.886006   \n\n                 dse  warning loo_scale  \nlog-normal  0.000000    False       log  \nexgaussian  2.742772    False       log  \ngamma       1.221702    False       log  \nnormal      2.213495    False       log  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>rank</th>\n      <th>loo</th>\n      <th>p_loo</th>\n      <th>d_loo</th>\n      <th>weight</th>\n      <th>se</th>\n      <th>dse</th>\n      <th>warning</th>\n      <th>loo_scale</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>log-normal</th>\n      <td>0</td>\n      <td>-612.674051</td>\n      <td>3.196749</td>\n      <td>0.000000</td>\n      <td>6.360410e-01</td>\n      <td>8.412337</td>\n      <td>0.000000</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>exgaussian</th>\n      <td>1</td>\n      <td>-613.321656</td>\n      <td>4.551529</td>\n      <td>0.647606</td>\n      <td>3.639590e-01</td>\n      <td>7.838731</td>\n      <td>2.742772</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>gamma</th>\n      <td>2</td>\n      <td>-614.527448</td>\n      <td>2.904233</td>\n      <td>1.853397</td>\n      <td>0.000000e+00</td>\n      <td>7.928021</td>\n      <td>1.221702</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>normal</th>\n      <td>3</td>\n      <td>-620.535490</td>\n      <td>3.511112</td>\n      <td>7.861439</td>\n      <td>1.287859e-14</td>\n      <td>8.886006</td>\n      <td>2.213495</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":45}],"execution_count":45},{"cell_type":"markdown","source":"- 第一列为索引，它列出了传递给 az.compare(.) 的模型名称。\n\n- rank 列：按照预测精度做的排名，值从0依次到模型总数，其中0代表最高精度。\n\n- loo 列：各模型 ELPD 值的列表，总是按照 ELPD 值从最好到最差排序。\n\n- p_loo 列：惩罚项的值列表，可以将其粗略地视为有效参数数量的估计值（但不要太认真）。此值可能低于具有更多结构的模型（如分层模型）中的实际参数数量，或者高于那些预测能力非常弱或严重错误指定的模型的实际参数数量。\n\n- d_loo 列：每个模型与排名第一的模型之间的 LOO 相对差。因此第一个模型始终取值为0  。\n\n- weight 列：分配给每个模型的权重。权重可以粗略地解释为在指定数据的条件下，是（参与比较的各模型中）该模型的概率。\n\n- se 列：ELPD 的标准误差。\n\n- dse 列：ELPD 相对差的标准误差。 dse 与 se 不一定相同，因为 ELPD 的不确定性在模型之间可能存在相关性。排名第一的模型 dse 值始终为0 。\n\n- warning 列：如果为True，表示这是一个警告，LOO 的近似估计不可靠。\n\n- loo_scale 列：估计值所用的尺度。默认为对数尺度。其他选项还包括：离差值尺度，即对数分值乘以-2，这会颠倒排序，较低的 ELPD 会更好；负对数尺度，即对数分值乘以-1，与离差值尺度一样，值越低越好。\n","metadata":{"id":"AF9DA5C23DFB4DFC9F3C5C440A280C3A","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"#### WAIC","metadata":{"id":"5DEC24C57D93492B888FBC3B4EE250D9","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"# 将4个模型的采样结果进行比较\ncompare_dict = {\"normal\": trace1, \"log-normal\": trace2, \"gamma\": trace3,\"exgaussian\":trace4}\n# 选择loo方法进行比较\ncomp = az.compare(compare_dict, ic='waic')\ncomp","metadata":{"id":"06E91DFD809B42C388C3D6DB2915D758","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n  \"The default method used to estimate the weights for each model,\"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n"},{"output_type":"execute_result","data":{"text/plain":"            rank        waic    p_waic    d_waic        weight        se  \\\nlog-normal     0 -612.655561  3.178259  0.000000  6.283489e-01  8.402979   \nexgaussian     1 -613.222991  4.452864  0.567430  3.716511e-01  7.809461   \ngamma          2 -614.520046  2.896831  1.864485  0.000000e+00  7.925004   \nnormal         3 -620.530952  3.506573  7.875390  6.772360e-15  8.883640   \n\n                 dse  warning waic_scale  \nlog-normal  0.000000     True        log  \nexgaussian  2.719900     True        log  \ngamma       1.218672     True        log  \nnormal      2.218972     True        log  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>rank</th>\n      <th>waic</th>\n      <th>p_waic</th>\n      <th>d_waic</th>\n      <th>weight</th>\n      <th>se</th>\n      <th>dse</th>\n      <th>warning</th>\n      <th>waic_scale</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>log-normal</th>\n      <td>0</td>\n      <td>-612.655561</td>\n      <td>3.178259</td>\n      <td>0.000000</td>\n      <td>6.283489e-01</td>\n      <td>8.402979</td>\n      <td>0.000000</td>\n      <td>True</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>exgaussian</th>\n      <td>1</td>\n      <td>-613.222991</td>\n      <td>4.452864</td>\n      <td>0.567430</td>\n      <td>3.716511e-01</td>\n      <td>7.809461</td>\n      <td>2.719900</td>\n      <td>True</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>gamma</th>\n      <td>2</td>\n      <td>-614.520046</td>\n      <td>2.896831</td>\n      <td>1.864485</td>\n      <td>0.000000e+00</td>\n      <td>7.925004</td>\n      <td>1.218672</td>\n      <td>True</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>normal</th>\n      <td>3</td>\n      <td>-620.530952</td>\n      <td>3.506573</td>\n      <td>7.875390</td>\n      <td>6.772360e-15</td>\n      <td>8.883640</td>\n      <td>2.218972</td>\n      <td>True</td>\n      <td>log</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":46}],"execution_count":46},{"cell_type":"markdown","source":"#### DIC & others","metadata":{"id":"AB229CA7F4CF47D094C426AC853CA97B","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"def comp_model(trace, model):\n    '''\n    trace: 未转化为inferencedata的采样结果\n    model: 模型\n    '''\n    dftrc_m = pm.trace_to_dataframe(trace, include_transformed=True) # 将trace对象转化为dataframe\n    trace = az.from_pymc3(trace) # 将trace对象转化为inferencedata对象\n    trc_logp = trace['log_likelihood']['y_obs'].to_dataframe().groupby(['chain','draw']).sum().reset_index()['y_obs'] #从模型中提取对数似然\n    \n    #dic\n    mean_deviance = -2 * trc_logp.mean(0) #计算所有参数对数似然的平均值\n    deviance_at_mean = -2 * model.logp(dftrc_m.mean(0).to_dict()) #计算参数平均值的对数似然\n    dic = 2 * mean_deviance - deviance_at_mean\n    \n    #bic\n    #deviance_at_mle = min(trc_logp) #对数似然最小的值\n    #parnum = 3 # 参数数量\n    #n = 3988 # 数据的样本数\n    #bic = -2 * deviance_at_mle +  2*parnum*np.log(n)\n    \n    #aic\n    #aic = -2 * deviance_at_mle +  2*parnum\n    \n    return dic\n    #,bic,aic","metadata":{"id":"F8ACDFA995524B1A8CF222A805F48957","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":47},{"cell_type":"code","source":"# 计算各模型dic的值\ndic2 = comp_model(trace=trace2_for_comp,model=LogNormal)\ndic4 = comp_model(trace=trace4_for_comp,model=exgaussian)\ndic3 = comp_model(trace=trace3_for_comp,model=Gamma)\ndic1 = comp_model(trace=trace1_for_comp,model=NormalModel)\ndics = {'gamma':dic3,'log-normal':dic2,'exgaussian':dic4,'normal':dic1}","metadata":{"id":"BF96BA27FE7F4A0189AF6990400D97F1","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n"}],"execution_count":48},{"cell_type":"code","source":"model_comp = pd.DataFrame({ \n                        'rank':az.compare(compare_dict, ic='waic')['rank'],\n                        'waic': az.compare(compare_dict, ic='waic')['waic'],\n                        'loo': az.compare(compare_dict, ic='loo')['loo'],\n                        })\nmodel_comp.loc['log-normal','dic'] =dics['log-normal']\nmodel_comp.loc['normal','dic'] =dics['normal']\nmodel_comp.loc['exgaussian','dic'] =dics['exgaussian']\nmodel_comp.loc['gamma','dic'] =dics['gamma']","metadata":{"id":"70CB60C7BEBA4A1C8B53015E0E395740","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":true,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n  \"The default method used to estimate the weights for each model,\"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n  \"The default method used to estimate the weights for each model,\"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:1407: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \nSee http://arxiv.org/abs/1507.04544 for details\n  \"For one or more samples the posterior variance of the log predictive \"\n/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n  \"The default method used to estimate the weights for each model,\"\n"}],"execution_count":49},{"cell_type":"code","source":"model_comp","metadata":{"id":"F0E59C78ACAB4CA79947388D98E4787A","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"            rank        waic         loo          dic\nlog-normal     0 -612.655561 -612.674051  1175.294384\nexgaussian     1 -613.222991 -613.321656  1194.861382\ngamma          2 -614.520046 -614.527448  1201.247458\nnormal         3 -620.530952 -620.535490  1212.743009","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>rank</th>\n      <th>waic</th>\n      <th>loo</th>\n      <th>dic</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>log-normal</th>\n      <td>0</td>\n      <td>-612.655561</td>\n      <td>-612.674051</td>\n      <td>1175.294384</td>\n    </tr>\n    <tr>\n      <th>exgaussian</th>\n      <td>1</td>\n      <td>-613.222991</td>\n      <td>-613.321656</td>\n      <td>1194.861382</td>\n    </tr>\n    <tr>\n      <th>gamma</th>\n      <td>2</td>\n      <td>-614.520046</td>\n      <td>-614.527448</td>\n      <td>1201.247458</td>\n    </tr>\n    <tr>\n      <th>normal</th>\n      <td>3</td>\n      <td>-620.530952</td>\n      <td>-620.535490</td>\n      <td>1212.743009</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":50}],"execution_count":50},{"cell_type":"markdown","source":"## Part 2: 广义线性模型(Generalized linear model, GLM)","metadata":{"id":"E39532601883414EBBF66B636BF81C9E","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"本节的目的在与：了解如何通过广义线性模型(Generalized linear model, GLM)拟合正确率等二元决策变量。\n\n重点在于：\n- 了解使用 Pymc 进行数据分析的完整 workflow\n- 了解因变量为正确率等二元决策变量(往往记录为0或1，1代表回答正确，0代表回答错误)的特征\n- 了解广义线性模型(Generalized linear model)中的伯努利(Bernoulli)分布和链接函数(link function)\n","metadata":{"id":"6B25238778D84118847CE06DE9047E8D","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"### 什么是GLM？","metadata":{"id":"13173C99D55F44AB8911996B1C0DF9EF","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"GLM 是一般线性模型的一种扩展。我们首先介绍一般线性模型的基础知识。\n\n在一般线性模型中，因变量被假定为服从正态分布 $y \\sim Normal(\\mu,sigma)$。\n- 其中，y为观测项；$\\mu$为预测项；sigma 为误差项。\n- 预测项展开为 $\\mu  = \\alpha + \\beta *x$。\n- 例子，比如用员工工龄(x)预测他们的工资(y)，其中x和y都为连续变量。\n\n\n![Image Name](https://cdn.kesci.com/upload/image/rll49b8jn9.png?imageView2/0/w/640/h/640)","metadata":{"id":"BF567BDD20DB4A4D82FB107EECF9F8D8","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"比如，我们这里男性编码为0，将女性编码为1。\n\n$Y = \\beta_0 + \\beta_1 X + \\epsilon$\n- 当虚拟变量赋值为X=0时， $E(Y) = \\beta_0$ 代表男性的平均家庭收入\n- 当虚拟变量赋值为X=1时，$E(Y) = \\beta_0+\\beta_1$ 代表女性的平均家庭收入\n- $\\beta_1$ 表示女性相对于男性的家庭收入的差值","metadata":{"id":"6D1C1DF77FFE43C78B561EA22ECC26F5","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"- 当自变量为二分变量时，t 检验是线性模型的特殊形式。包括独立样本t检验和配对样本t检验。\n- 当自变量为超过俩水平的分类变量时，方差分析是线性模型的特殊形式。包括方差分析和重复测量方差分析。\n\n例子中的性别为二分变量，对应独立样本t检验。\n\n与独立样本t检验不同，在线性模型中检查两组之间的差异相当于检查回归系数 $b_1 $的显著性。\n\n![Image Name](https://cdn.kesci.com/upload/image/rloa41zjxa.png?imageView2/0/w/640/h/640)\n","metadata":{"id":"C3AAC8E60AE54BE198EB83BC70FAF859","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"**当因变量为离散变量**\n\n比如，因变量为答题正确率，其中1代表回答正确，0代表回答错误。\n\n正确率为离散变量，并不服从正态分布，而是服从伯努利(Bernoulli)分布。\n\n\n![Image Name](https://cdn.kesci.com/upload/image/rloa62fn8w.png?imageView2/0/w/960/h/960)\n\n\n![](https://docs.pymc.io/en/v4.3.0/_images/pymc-Bernoulli-1.png)","metadata":{"id":"D80B8735D4AA445B8ECD14F8CECCD737","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"广义线性模型(Generalized linear model，GLM)的特点：\n\n| 一般线性模型 | 广义线性模型 | \n|---|---|\n| $y \\sim Normal(\\mu,sigma)$ | $y \\sim dist(p)$ |\n| $\\mu = \\alpha + \\beta *x$ | $p = g(\\mu)$|\n\n\n- 首先，GLM 可以将 $y \\sim Normal(\\mu,sigma)$ 扩展为 $y \\sim Bernoulli(p)$ ，使得因变量y服从伯努利分布。\n- 同样，参数 $p$ 可以与自变量联系在一起， $p  = \\alpha + \\beta * x$。\n- 需要注意的是，由于$p$ 的取值范围在[0, 1]，而 $\\alpha + \\beta * x$ 的范围为 $(-\\infty, +\\infty)$。我们需要通过**链接函数g()** 将 $\\alpha + \\beta * x$  映射到 $p$ 所在的范围。\n\n\n","metadata":{"id":"DB1B352289A041ACB3C0E05E82394D98","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"链接函数的具体转化过程，以逻辑(logit)回归为例：\n1. 令 $z = \\alpha + \\beta *x$，$\\mu$的范围为 $(-\\infty, +\\infty)$。\n2. $p = g(z)$，其中 g() 为链接函数，输出结果 p 的范围为 $(0,1)$。\n3.  最后将 p 输入到分布函数中 $y \\sim Bernoulli(p)$。\n\n\n![Image Name](https://cdn.kesci.com/upload/image/rloa6zyf5a.png?imageView2/0/w/600/h/600)","metadata":{"id":"F77E83AC8951422E84D349179CD09B76","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"### Workflow\n\n在了解GLM的基础知识后，我们通过实际的例子来体验PyMC建模的完整workflow。\n\n![Image Name](https://cdn.kesci.com/upload/image/rkvikqg9q6.png?imageView2/0/w/650/h/650)","metadata":{"id":"1C7624EF3B504B35AFED28D059769B4D","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"### (1) 提出研究问题\n\n假如研究问题为：不同刺激条件下**正确率**的差异。\n\n![](https://www.researchgate.net/profile/Ata_Akin/publication/281167153/figure/download/fig1/AS:391418049777669@1470332743733/Three-different-stimulus-conditions-in-the-Stroop-task-neutral-congruent-and.png?_sg=ibeklp8QZ2sbyR29ZZxbOgfS--_RjcKP_uVY36qBahzEJlnMLYPxQyzgYT2Au85eDBClhLqol0A)\n\n\n\n图片来源：https://www.researchgate.net/publication/281167153_Similarity_analysis_of_functional_connectivity_with_functional_near-infrared_spectroscopy/figures?lo=1&utm_source=bing&utm_medium=organic","metadata":{"id":"6ADD8C412C3742788AC1504A7DEA4ED2","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"### (2) 数据收集","metadata":{"id":"A607902A947049F997E399A6C8631382","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"#加载需要使用的库\n%matplotlib inline\nimport numpy as np \nfrom scipy import stats\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport arviz as az\nimport pymc3 as pm\n# 随机数种子，确保随后生成的随机数相同\nnp.random.seed(123)  ","metadata":{"id":"54EC70848A8E4DED9FAAD4D46A2F8905","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":51},{"cell_type":"markdown","source":"这里我们使用的数据来自 Eisenberg et al (2019)。\n\n为了简化问题，我们仅考虑有一个被试的数据。其中：\n- worker_id 为被试编号。\n- correct 为被试在 stroop 任务中每个试次判断的正确性，其中1代表判断正确，0代表判断错误。\n- condition 为刺激的类别，congruent为颜色和字意一致，incongruent为颜色和字意不一致。","metadata":{"id":"42472D58A17C4A34811A6E4397459D0F","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"# 加载数据\ndata = pd.read_csv(\"/home/mw/input/data4676/stroop.csv\")\n# 选取第一个被试在正式实验(test)中的数据\ndata = data[(data.worker_id == \"s001\") & (data.exp_stage == \"test\")]\n# 选取数据中的判断正确率，刺激条件，和被试编号\ndata = data[[\"worker_id\",\"correct\",\"condition\"]]\n# 重置每个试次的编号\ndata.reset_index(inplace=True,drop=True)","metadata":{"id":"D73129F702E14F9A826935B969C52FA9","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":52},{"cell_type":"code","source":"data.head()","metadata":{"id":"86126A55F30F43F9813C6F48CD60C1CB","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"  worker_id  correct  condition\n0      s001      1.0  congruent\n1      s001      0.0  congruent\n2      s001      1.0  congruent\n3      s001      1.0  congruent\n4      s001      1.0  congruent","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>worker_id</th>\n      <th>correct</th>\n      <th>condition</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>s001</td>\n      <td>1.0</td>\n      <td>congruent</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>s001</td>\n      <td>0.0</td>\n      <td>congruent</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>s001</td>\n      <td>1.0</td>\n      <td>congruent</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>s001</td>\n      <td>1.0</td>\n      <td>congruent</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>s001</td>\n      <td>1.0</td>\n      <td>congruent</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":53}],"execution_count":53},{"cell_type":"code","source":"data.groupby('condition').correct.describe() ","metadata":{"id":"477B65276EE74200867E1AF1BFF2211F","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"             count    mean       std  min  25%  50%  75%  max\ncondition                                                    \ncongruent     48.0  0.8750  0.334219  0.0  1.0  1.0  1.0  1.0\nincongruent   48.0  0.8125  0.394443  0.0  1.0  1.0  1.0  1.0","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>count</th>\n      <th>mean</th>\n      <th>std</th>\n      <th>min</th>\n      <th>25%</th>\n      <th>50%</th>\n      <th>75%</th>\n      <th>max</th>\n    </tr>\n    <tr>\n      <th>condition</th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>congruent</th>\n      <td>48.0</td>\n      <td>0.8750</td>\n      <td>0.334219</td>\n      <td>0.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>incongruent</th>\n      <td>48.0</td>\n      <td>0.8125</td>\n      <td>0.394443</td>\n      <td>0.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":54}],"execution_count":54},{"cell_type":"code","source":"data.groupby(['condition']).correct.mean().plot.bar()\nplt.show()","metadata":{"id":"3401F4AD272944C78FE42B793CFD6539","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/3401F4AD272944C78FE42B793CFD6539/rlox10fkhe.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":55},{"cell_type":"markdown","source":"### (3) 选择模型","metadata":{"id":"A9889876AADA4811B1565317B8ACAFC0","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"在我们的例子中，由于因变量(反应的正确性)不是连续变量，因此预测变量不服从正态分布。\n\n二分变量通常采用伯努利(Bernoulli)分布进行分析，因此我们需要广义线性模型(Generalized linear model，GLM)来扩展一般线性模型：\n- 首先，GLM 可以将 $y \\sim Normal(\\mu,sigma)$ 扩展为 **$y \\sim Bernoulli(p)$** ，使得因变量y服从伯努利分布。\n- 同样，参数 $p$ 可以与自变量联系在一起， $p  = \\alpha + \\beta * x$。\n- 需要注意的是，$p$ 的取值范围在[0, 1]，而 $\\alpha + \\beta * x$ 的范围为 $(-\\infty, +\\infty)$。我们需要通过**链接函数** 将 $\\alpha + \\beta * x$  映射到 $p$ 所在的范围。\n\t1. 令 $z = \\alpha + \\beta *x$，$\\mu$的范围为 $(-\\infty, +\\infty)$。\n\t2. $p = g(z)$，其中 g() 为链接函数，输出结果 $p$ 的范围为 $(0,1)$。\n\t3.  最后将 $p$ 输入到分布函数中 $y \\sim Bernoulli(p)$。","metadata":{"id":"7A70297376F74E388F4F7D1FC92E18D7","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"### (4)选择先验","metadata":{"id":"FC76212253944A3C845433F6F2B699F0","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"# 将‘condition’进行编码，其中一致条件(congruent)编码为0，不一致条件(incongruent)编码为1。\ndata.condition = data.condition.map({'incongruent':1,'congruent':0})","metadata":{"id":"17D947CB6C4D43BB8C603A91D8F7C316","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":56},{"cell_type":"code","source":"# 在pymc3中，pm.Model()定义了一个新的模型对象，这个对象是模型中随机变量的容器\n# 在python中，容器是一种数据结构，是用来管理特殊数据的对象\n# with语句定义了一个上下文管理器，以 linear_model 作为上下文，在这个上下文中定义的变量都被添加到这个模型\nwith pm.Model() as GLM_model:\n    # 设定先验分布: \n    alpha = pm.Normal('alpha',mu=0,sd=1)\n    beta = pm.Normal('beta',mu=0,sd=1)\n    # x为自变量，是之前已经载入的数据\n    x = pm.Data(\"x\", data['condition'])\n    # 通过链接函数对参数进行转换\n    z = alpha + beta * x                            # 对应步骤1\n    p = pm.Deterministic(\"p\", pm.math.invlogit(z))  # 对应步骤2\n\n    # 先验预测检查\n    prior_checks = pm.sample_prior_predictive(samples=50)","metadata":{"id":"6BCEBA5F827348BA826278FF87AE2D46","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":57},{"cell_type":"code","source":"az.plot_density(\n    {'alpha':prior_checks['alpha'],\n    'beta':prior_checks['beta']}\n    )\nplt.show()","metadata":{"id":"0A87F33FAA004F09B298629449E61FFE","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 993.6x331.2 with 2 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/0A87F33FAA004F09B298629449E61FFE/rlox10zt5m.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":58},{"cell_type":"code","source":"az.plot_density(\n    {'p':prior_checks['p']}\n    )\nplt.show()","metadata":{"id":"F327DDA4CA9A4D6083A80F16B56A3D66","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/F327DDA4CA9A4D6083A80F16B56A3D66/rlox10tgfh.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":59},{"cell_type":"markdown","source":"结果发现，通过**链接函数**转换后的p值范围为 0到1。","metadata":{"id":"681CAA1FEDE64E5497F77C5EFD168F6B","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"### (5) 拟合数据","metadata":{"id":"6EF2BFFA7F454291A705B45D91343399","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"首先定义 GLM 模型：\n- 其中 alpha 和 beta 为模型参数，$\\alpha + \\beta * x$。\n- x 为自变量刺激条件(condition), 0代表一致条件，1代表不一致条件。\n- 通过链接函数对参数进行转换。\n    - 首先令 $\\mu = \\alpha + \\beta * x$\n    - 然后通过链接函数 pm.math.invlogit(mu)，计算出 $p$。\n    注意，这里我们选择 logit 的反函数 invlogit作为链接函数。该链接函数使得 $p$ 的范围为 $(0,1)$。\n-  最后将 $p$ 输入到分布函数中 $y \\sim Bernoulli(p)$。","metadata":{"id":"D5F0A7F2C63D4D4E8415C468F15BF07E","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with pm.Model() as GLM_model:\n    # 在pymc3中，pm.Model()定义了一个新的模型对象，这个对象是模型中随机变量的容器\n    # 在python中，容器是一种数据结构，是用来管理特殊数据的对象\n    # with语句定义了一个上下文管理器，以 linear_model 作为上下文，在这个上下文中定义的变量都被添加到这个模型\n    # 定义先验\n    alpha = pm.Normal('alpha', mu=0, sd=1)\n    beta = pm.Normal('beta', mu=0, sd=1)\n    # x为自变量，是之前已经载入的数据\n    x = pm.Data(\"x\", data['condition'])\n    # 线性模型：mu是确定性随机变量，这个变量的值完全由右端值确定\n    z = alpha + beta * x                            # 对应步骤1\n    p = pm.Deterministic(\"p\", pm.math.invlogit(z))  # 对应步骤2\n    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    y_obs = pm.Bernoulli(\"y_obs\", p=p, observed=data[\"correct\"])  # 对应步骤3","metadata":{"id":"E5BD75E7001346C394E1C251BDA49425","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":60},{"cell_type":"code","source":"# 展示模型结构\npm.model_to_graphviz(GLM_model)","metadata":{"id":"12DF7FF82BF044C5817107FC51005663","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7f2f621a9710>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/12DF7FF82BF044C5817107FC51005663/rlox1125zc.svg\">"},"metadata":{},"execution_count":61}],"execution_count":61},{"cell_type":"markdown","source":"注意：由于应用 GLM 模型时往往都会使用到链接函数，为了减轻使用者的工作量，在 pymc中可以通过设定 将 `pm.Bernoulli(\"y_obs\", p=p)` 的设定改为 `pm.Bernoulli(\"y_obs\", logit_p=p)` 。 完整代码如下：","metadata":{"id":"21E4758E9B7347D388E2012B438682E6","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with pm.Model() as GLM_model:\n    # 定义先验\n    alpha = pm.Normal('alpha', mu=0, sd=1)\n    beta = pm.Normal('beta', mu=0, sd=1, shape=1)\n    # x为自变量，是之前已经载入的数据\n    x = pm.Data(\"x\", data['condition'])\n    # 线性模型：mu是确定性随机变量，这个变量的值完全由右端值确定\n    p = pm.Deterministic(\"p\", alpha + beta * x)\n    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data[\"correct\"])","metadata":{"id":"A1DBBED0269842A4816495F01387C9B5","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":62},{"cell_type":"markdown","source":"### (6)采样过程诊断\n\n如果使用MCMC对后验进行近似，则需要首先对MCMC过程进行评估。\n\n* 是否收敛；\n* 是否接近真实的后验。\n\n对采样过程的评估我们会采用目视检查或rhat这个指标","metadata":{"id":"9F3B001B8A8A4062B623829660B1F32E","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with GLM_model :\n    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n    trace = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)","metadata":{"id":"3BD7E772AF6B4FDD9510B506DB02E7BF","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [beta, alpha]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 5 seconds.\n"}],"execution_count":63},{"cell_type":"code","source":"az.plot_trace(trace, var_names=['alpha','beta'])\nplt.show()","metadata":{"id":"731C1B5025A44B03851429397F824611","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 864x288 with 4 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/731C1B5025A44B03851429397F824611/rlox188dnp.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":64},{"cell_type":"code","source":"az.summary(trace, var_names=['alpha','beta'], kind=\"diagnostics\")","metadata":{"id":"F1C4D1787A2549AF880E775EBF9EE51D","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"         mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat\nalpha        0.010    0.007    1314.0    1578.0    1.0\nbeta[0]      0.013    0.010    1267.0    1500.0    1.0","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>alpha</th>\n      <td>0.010</td>\n      <td>0.007</td>\n      <td>1314.0</td>\n      <td>1578.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>beta[0]</th>\n      <td>0.013</td>\n      <td>0.010</td>\n      <td>1267.0</td>\n      <td>1500.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":65}],"execution_count":65},{"cell_type":"markdown","source":"### (7)模型诊断\n\n在MCMC有效的前提下，需要继续检验模型是否能够较好地拟合数据。\n\n我们会使用后验预测分布通过我们得到的参数生成一批模拟数据，并将其与真实数据进行对比。","metadata":{"id":"48E2254DF40E43238B0DCDBB19BADD78","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"# 后验预测分布的计算仍在容器中进行\nwith GLM_model:\n    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n    ppc_y = pm.sample_posterior_predictive(trace.posterior) \n#将ppc_y转化为InferenceData对象合并到trace中\naz.concat(trace, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)","metadata":{"id":"C79C78938F80443A9B700DB7792CDA63","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"update_display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n    <div>\n      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      100.00% [4000/4000 00:03&lt;00:00]\n    </div>\n    "},"metadata":{}},{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n"}],"execution_count":66},{"cell_type":"code","source":"# 绘制后验预测分布\naz.plot_ppc(trace)\nplt.show()","metadata":{"id":"23CCDDD003D244A48135386EFAB8F7F0","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/IPython/core/pylabtools.py:151: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n  fig.canvas.print_figure(bytes_io, **kw)\n"},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/23CCDDD003D244A48135386EFAB8F7F0/rlox1n5ygd.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":67},{"cell_type":"markdown","source":"### (8)模型比较","metadata":{"id":"228CA1391AF84CE187FBFFF3BDACEF63","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"当采样诊断与模型诊断说明模型是否可用后，我们可以通过模型检验来验证我们的研究问题：在不同刺激条件下(一致 vs. 不一致)个体正确率是否存在差异。\n\n![Image Name](https://cdn.kesci.com/upload/image/rkm3pw954u.png?imageView2/0/w/960/h/960)","metadata":{"id":"D5CEDE077C1D48628B59F5D67D9C1F26","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"我们可以定义一个不考虑自变量影响的模型 `GLM_null_model`。如果之前的模型拟合优度好于该模型，那么说明自变量对模型存在影响。","metadata":{"id":"236F899260DC460B8F957A5CA47C2058","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"with pm.Model() as GLM_null_model:\n    # 定义先验\n    p = pm.Uniform('p', 0, 1)  # 由于没有考虑自变量的影响，因此我们可以直接假设参数p服从0到1的均匀分布。\n    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    y_obs = pm.Bernoulli('y_obs', p=p, observed=data['correct'])\n\n    trace2 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)","metadata":{"id":"E93A626AEBC64B1F8D36F83B09CDA081","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [p]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 3 seconds.\n"}],"execution_count":68},{"cell_type":"code","source":"############################\n# 练习\n# 要求：完成对 GLM_null_model 的采样过程诊断与模型诊断。\n############################\n\n# 绘制各参数的采样情况\n# Tips: 使用 az.plot_trace() 函数可以绘制 trace 图；使用 az.summary() 可以得到诊断统计结果\n\n\n# 模型诊断\n# Tips: 使用 pm.sample_posterior_predictive 可以进行后验预测检验；使用 az.plot_ppc() 可以得到后验预测检验图\n","metadata":{"id":"196F0C31636C4BA9874E0B114592537A","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[],"execution_count":69},{"cell_type":"markdown","source":"当对 GLM_null_model 进行同样的检验，我们可以正式进行模型比较了。","metadata":{"id":"DC9631045212432F810EB1AA05FF7F5A","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"# 将三个模型的采样结果进行比较\ncompare_dict = {\"GLM_model\": trace, \"GLM_null_model\": trace2}\n# 选择loo方法进行比较\ncomp = az.compare(compare_dict, ic='loo')\ncomp","metadata":{"id":"D7E2770EE0F046878DB9DAB31507489E","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n  \"The default method used to estimate the weights for each model,\"\n"},{"output_type":"execute_result","data":{"text/plain":"                rank        loo     p_loo     d_loo  weight        se  \\\nGLM_null_model     0 -42.584941  0.954592  0.000000     1.0  5.963029   \nGLM_model          1 -43.136095  1.660994  0.551154     0.0  5.839979   \n\n                     dse  warning loo_scale  \nGLM_null_model  0.000000    False       log  \nGLM_model       0.371372    False       log  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>rank</th>\n      <th>loo</th>\n      <th>p_loo</th>\n      <th>d_loo</th>\n      <th>weight</th>\n      <th>se</th>\n      <th>dse</th>\n      <th>warning</th>\n      <th>loo_scale</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>GLM_null_model</th>\n      <td>0</td>\n      <td>-42.584941</td>\n      <td>0.954592</td>\n      <td>0.000000</td>\n      <td>1.0</td>\n      <td>5.963029</td>\n      <td>0.000000</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>GLM_model</th>\n      <td>1</td>\n      <td>-43.136095</td>\n      <td>1.660994</td>\n      <td>0.551154</td>\n      <td>0.0</td>\n      <td>5.839979</td>\n      <td>0.371372</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":70}],"execution_count":70},{"cell_type":"markdown","source":"结果显示，`GLM_null_model` 模型的拟合度好于 `GLM_model`，说明不存在充分的证据表明不同刺激条件会影响个体判断的正确率。","metadata":{"id":"7C007AE6B6E9435781F86808C0FB5997","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"### (9)统计推断","metadata":{"id":"A6E2192CE2FA4E438F867CCE1EDB9726","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"我们可以进一步通过统计推断印证模型比较的结果。","metadata":{"id":"C52ABCE8AA3340F19E4C6DF68952C4C7","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"az.plot_posterior(trace, var_names=['beta'])\nplt.show()","metadata":{"id":"A8A1EEC2E1C54598855026A223ACEDBF","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/A8A1EEC2E1C54598855026A223ACEDBF/rlox1tjmfq.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":71},{"cell_type":"markdown","source":"参数 $\\beta$ 反应了两个条件下正确率的差异。我们可以看到，该参数的后验分布的大部分包括0，说明支持两个条件下正确率存在差异的证据不足。","metadata":{"id":"B0193DCF428D462E82BAF2D0E3DEBD38","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"我们进一步查看两个参数的情况：","metadata":{"id":"DA6F68ADBF7C40F59D0DEFCF39D3FA52","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"az.summary(trace, var_names=['alpha','beta'])","metadata":{"id":"FD7DCADA005F4EC49CA9C78D199D9337","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"          mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\nalpha    1.706  0.361   1.047    2.374      0.010    0.007    1314.0   \nbeta[0] -0.193  0.477  -1.103    0.716      0.013    0.010    1267.0   \n\n         ess_tail  r_hat  \nalpha      1578.0    1.0  \nbeta[0]    1500.0    1.0  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>alpha</th>\n      <td>1.706</td>\n      <td>0.361</td>\n      <td>1.047</td>\n      <td>2.374</td>\n      <td>0.010</td>\n      <td>0.007</td>\n      <td>1314.0</td>\n      <td>1578.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>beta[0]</th>\n      <td>-0.193</td>\n      <td>0.477</td>\n      <td>-1.103</td>\n      <td>0.716</td>\n      <td>0.013</td>\n      <td>0.010</td>\n      <td>1267.0</td>\n      <td>1500.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":72}],"execution_count":72},{"cell_type":"markdown","source":"结果发现，两个参数值的范围不太“正常”。\n\n还记得 **链接函数**吗？\n- 链接函数 g() 将 $\\alpha + \\beta * x$ 的范围从 $(-\\infty, +\\infty)$ 转换为 (0,1)\n- 同时，其中的参数 $\\alpha$ 和 $\\beta$也被转换了，只不过他们从 (0,1) 转换为  $(-\\infty, +\\infty)$，所以  $\\alpha$ 大于1，并且 $\\beta$小于0。\n- 为了他们转换回来，我们需要使用 logit 函数， $p = \\frac{1}{1+e^θ}$。\n具体代码如下：","metadata":{"id":"7E9B9DC8B5CB481E89D2167099E3A3AA","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"p_congruent = 1 / (1 + np.exp(-trace.posterior[\"alpha\"].mean())).to_pandas()\np_incongruent = 1 / (1 + np.exp(-(trace.posterior[\"beta\"].mean()+trace.posterior[\"alpha\"].mean()))).to_pandas()\nprint(\"alpha(一致条件) = \",p_congruent, \"\\n alpha+beta(不一致条件) = \", p_incongruent)","metadata":{"id":"72062F4C925E425789BC962D35C47709","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stdout","text":"alpha(一致条件) =  0.8462798776599924 \n alpha+beta(不一致条件) =  0.8194362023007759\n"}],"execution_count":73},{"cell_type":"markdown","source":"转换后可以发现，虽然一致条件的正确率略高于不一致条件，但这个差异并不具有统计学意义。","metadata":{"id":"ECA366CCCD4449D289E2EA2C7805BD4D","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"markdown","source":"值得注意的是，模型预测的正确率与实际数据的正确率存在差异。\n- 对于一致条件刺激，模型预测值为0.842，**低于**实际正确率为0.875。\n- 对于不一致条件刺激，模型预测值为0.823，**高于**实际正确率为0.813。\n- 模型预测两个条件的正确的差为(一致条件-不一致条件) = 0.02，而真实正确率的差异 = 0.06\n\n可见，模型预测的效应更小。这是因为我们设置先验时，认为条件间的差异(即beta参数)服从均值为0，标准差为1的正态分布。","metadata":{"id":"392B5AF334684B9581977F9DD24D9DE1","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"data.groupby('condition').correct.describe() ","metadata":{"id":"07AF647FA2FF49368D7C5DEC0B69B969","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"           count    mean       std  min  25%  50%  75%  max\ncondition                                                  \n0           48.0  0.8750  0.334219  0.0  1.0  1.0  1.0  1.0\n1           48.0  0.8125  0.394443  0.0  1.0  1.0  1.0  1.0","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>count</th>\n      <th>mean</th>\n      <th>std</th>\n      <th>min</th>\n      <th>25%</th>\n      <th>50%</th>\n      <th>75%</th>\n      <th>max</th>\n    </tr>\n    <tr>\n      <th>condition</th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>48.0</td>\n      <td>0.8750</td>\n      <td>0.334219</td>\n      <td>0.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>48.0</td>\n      <td>0.8125</td>\n      <td>0.394443</td>\n      <td>0.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":74}],"execution_count":74},{"cell_type":"markdown","source":"假如我们存在先验知识，知道不一致条件的正确率低于一致条件，而不是他们可能相等。\n\n我们可以设置一个有信息的先验，比如假设这个差异为 0.6，即 beta = 0.6。\n\n需要注意的是，由于链接函数的存在，我们需要把 beta进行 logit转化，得到转化后的 beta 为 0.41。 \n\n因此，我们假定beta的先验分布服从 均值为 0.4，标准差为0.2的正态分布。使得beta大部分值都大于0 (相对于之前得到的beta = -0.134 < 0)。","metadata":{"id":"5F52E25DBD524D69876E7E0FAB173524","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}},{"cell_type":"code","source":"p = 0.6\nbeta = np.log(p/(1-p))\nprint(\"转化后的beta:\", beta)","metadata":{"id":"06F1D8FC005048619239D390C4689E5F","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stdout","text":"转化后的beta: 0.4054651081081642\n"}],"execution_count":75},{"cell_type":"code","source":"with pm.Model() as GLM_model2:\n    # 定义先验\n    alpha = pm.Normal('alpha',mu=0,sd=1)\n    beta = pm.Normal('beta',mu=0.4,sd=0.2) # 我们假定beta的先验分布服从 均值为 0.4，标准差为0.2的正态分布\n    # x为自变量，是之前已经载入的数据\n    x = pm.Data(\"x\", data['condition'])\n    # 线性模型：mu是确定性随机变量，这个变量的值完全由右端值确定\n    z = alpha - beta * x                  # 这里使用减号是因为一致条件比不一致条件的正确率高          \n    p = pm.Deterministic(\"p\", pm.math.invlogit(z))  \n    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    y_obs = pm.Bernoulli(\"y_obs\", p=p, observed=data[\"correct\"])  \n    trace3 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)","metadata":{"id":"090C4799FB654E6C850007C85F2EA986","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stderr","text":"Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [beta, alpha]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 5 seconds.\n"}],"execution_count":76},{"cell_type":"code","source":"p_congruent = 1 / (1 + np.exp(-trace3.posterior[\"alpha\"].mean())).to_pandas()\np_incongruent = 1 / (1 + np.exp(-(trace3.posterior[\"beta\"].mean()+trace3.posterior[\"alpha\"].mean()))).to_pandas()\nprint(\"alpha(一致条件) = \", p_congruent, \"\\n alpha+beta(不一致条件) = \", p_incongruent)","metadata":{"id":"A63FD397ED2A439E830E39926B15C2E7","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true},"outputs":[{"output_type":"stream","name":"stdout","text":"alpha(一致条件) =  0.855171478001516 \n alpha+beta(不一致条件) =  0.8956504809430272\n"}],"execution_count":77},{"cell_type":"markdown","source":"- 上一个模型模型预测两个条件的正确的差(一致条件-不一致条件) = 0.02，而真实正确率的差异 = 0.06。\n- 该模型预测两个条件的正确的差 = 0.05。已经非常接近真实差异。\n\n该结果说明了先验设置对于结果的影响。\n\n感兴趣的同学可以在课后尝试不同的先验设置，并进行相应的的模型诊断、模型比较和统计推断的练习。","metadata":{"id":"E1AA58B827FD473EACB16C40B0697A9B","notebookId":"637b29c276064fb3fe3e01b5","jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true}}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python","nbconvert_exporter":"python","file_extension":".py","version":"3.5.2","pygments_lexer":"ipython3"}},"nbformat":4,"nbformat_minor":0}